%%
%Ekman transport
%With amur,January-mean
%cross component vertical integration; integration within the larger Kv area
clc;clear;

%%
%read data
path(path,'/Users/yuannan/Desktop/sea ice in the okhotsk/in japan/gsw_matlab_v3_06_11');
name_file={'noriver','without-amur','0.25amur','0.5amur','with-amur','2amur','withamuronly',...
           'with-amur-5-ave','noriver-5-ave','without-amur-5-ave','withamuronly-5-ave','tide_JAN_monthly_mean'};
name_month={'Oct','Nov.','Dec.','Jan.','Feb.'}; 
file_seek=[9 10,11,0,1];
num_loop=1;
r1=[130:300];r2=[305:360];
name_var={'ahv','u','v','s','t','amv','wzc'};
name_var1={'ahv1','u1','v1','s1','t1','amv1','wzc1'};
name_var2={'ahv2','u2','v2','s2','t2','amv2','wzc2'};

cd /Volumes/Samsung_T5/program
read_parameter

for num_month=4
    for num_name_file=10

eval(['cd /Volumes/Samsung_T5/results_',name_file{num_name_file},'/data']);

%read data
for i=1:7
    eval(['fid',num2str(i),'=fopen(','''',name_var{i},'.grd','''',',','''','r','''',');']);
    if num_name_file==2 || num_name_file==5 || num_name_file==8
        eval(['fseek(fid',num2str(i),',800*800*84*4*','file_seek(1,num_month)',',-1);']); %for no fresh
    elseif num_name_file==1 || num_name_file==3 || num_name_file==4 || num_name_file==6 || num_name_file==7
        %eval(['fseek(fid',num2str(i),',800*800*84*4*12*4+800*800*84*4*','file_seek(1,num_month)',',-1);']); %for no fresh 查第5年数据
       eval(['fseek(fid',num2str(i),',800*800*84*4*12*9+800*800*84*4*','file_seek(1,num_month)',',-1);']); %for no fresh 查第10年数据
    end
    %eval(['fseek(fid',num2str(i),',800*800*84*4*12*9+800*800*84*4*','file_seek',',-1);']); %for no fresh
    eval([name_var{i},'=fread(fid',num2str(i),',','800*800*84',',','''','real*4','''',');']);
    eval([name_var1{i},'=reshape(',name_var{i},',800,800,84);']);
    eval([name_var1{i},'(',name_var1{i},'<-10000)=nan;']);
    eval(['fclose(fid',num2str(i),');']);
    eval([name_var2{i},'=',name_var1{i},'(r1,r2,:);']);
end

p(1:171,1:56,1:84)=10.1325;
rho2=gsw_rho(s2,t2,p)-1000;

u3=permute(u2,[2,1,3]);
v3=permute(v2,[2,1,3]);
rho3=permute(rho2,[2,1,3]);
amv3=permute(amv2,[2,1,3]);

% for z=1:84
%     u4(:,:,z)=u3(:,:,85-z);
%     v4(:,:,z)=v3(:,:,85-z);
%     rho4(:,:,z)=rho3(:,:,85-z);
%     amv4(:,:,z)=amv3(:,:,85-z);
% end
%%for last 5 years mean
for z=1:84
    u4(:,:,z)=u3(:,:,z);
    v4(:,:,z)=v3(:,:,z);
    rho4(:,:,z)=rho3(:,:,z);
    amv4(:,:,z)=amv3(:,:,z);
end

%partial rho/partial x and partial rho/partial y(gradient of rho)
for z=1:84
    %   [rho_x(:,:,z),rho_y(:,:,z)]=gradient(rho4(:,:,z),x1_dis,y1_dis);
       [rho_x(:,:,z),rho_y(:,:,z)]=gradient(rho4(:,:,z));
end
for z=1:84
    rho_x(:,:,z)=rho_x(:,:,z)./x_dis;
    rho_y(:,:,z)=rho_y(:,:,z)./y_dis;
end
%%
%cross- and along- direction
%coastal thermal wind
dens_sum_1(1:56,1:171)=0;
dens_sum_2(1:56,1:171)=0;
num_dens(1:56,1:171)=1;
dens_ave_1(1:56,1:171)=0;

for i=1:56
    for j=1:171
            for kk=2:84
                if isnan(rho4(i,j,kk))
                    continue
                else
                    dens_sum_2(i,j)=(rho4(i,j,kk)+rho4(i,j,kk-1))./2.0.*(depth(1,kk)-depth(1,kk-1));
                    dens_sum_1(i,j)=dens_sum_1(i,j)+dens_sum_2(i,j);
                    num_dens(i,j)=num_dens(i,j)+1;
                end
            end
        dens_ave_1(i,j)=dens_sum_1(i,j)./depth(1,num_dens(i,j));
    end
end

[densgx,densgy]=gradient(dens_ave_1);
densgx=densgx./x_dis;densgy=densgy./y_dis;
length_n=sqrt(densgx.*densgx+densgy.*densgy);
densgx=densgx./length_n;densgy=densgy./length_n; %gradient direction small--->large; right%downward:positive gradient
densay=sqrt(1./((densgy.*densgy)./(densgx.*densgx)+1));%tangent direction
densax=(-1.*densay.*densgy)./densgx;
%latitude 53.5N
for j=1:num_size1    
    gx_lat(j,1)=gx_new(test3(j,1),test4(j,1));
    gy_lat(j,1)=gy_new(test3(j,1),test4(j,1));  
    ax_lat(j,1)=ax_new(test3(j,1),test4(j,1));  
    ay_lat(j,1)=ay_new(test3(j,1),test4(j,1)); 
    densgx_lat(j,1)=densgx(test3(j,1),test4(j,1));  
    densgy_lat(j,1)=densgy(test3(j,1),test4(j,1));  
    densax_lat(j,1)=densax(test3(j,1),test4(j,1));  
    densay_lat(j,1)=densay(test3(j,1),test4(j,1));     
end
% scale_1=[39:43];%53.5N
% gx_lat(scale_1,1)=densgx_lat(scale_1,1);
% gy_lat(scale_1,1)=densgy_lat(scale_1,1);
% ax_lat(scale_1,1)=densax_lat(scale_1,1);
% ay_lat(scale_1,1)=densay_lat(scale_1,1);

%%
%geostrophic current
%latitude 53.5N
for j=1:num_size1
    amv_lat(j,:)=amv4(test3(j,1),test4(j,1),:);
    u_lat(j,:)=u4(test3(j,1),test4(j,1),:);
    v_lat(j,:)=v4(test3(j,1),test4(j,1),:);
    taupx(j,:)=-1*rho_y(test3(j,1),test4(j,1),:);
    taupy(j,:)=rho_x(test3(j,1),test4(j,1),:);
end
amv_lat=amv_lat.*0.0001;
u_lat=u_lat.*0.01;
v_lat=v_lat.*0.01;
speed=sqrt(u_lat.*u_lat+v_lat.*v_lat);
rho_constant=1023.6; %kg/m3
f_col=10e-5;%s^(-1)
g=9.8; %m/s^2
%Ekman thickness
deta_e=sqrt(amv_lat./(f_col/2));


%geostrophic current vertical shear(interior)
for z=1:84
    along(:,z)=u_lat(:,z).*ax_lat(:,1)+v_lat(:,z).*ay_lat(:,1);
    cross(:,z)=u_lat(:,z).*gx_lat(:,1)+v_lat(:,z).*gy_lat(:,1);
    taup_a(:,z)=-1*(taupx(:,z).*ax_lat(:,1)+taupy(:,z).*ay_lat(:,1)).*deta_e(:,z).*deta_e(:,z).*g./2; %cross- geostrophic stress
end

for j=1:num_size1
    for z=2:84
        along_z(j,z-1)=(along(j,z-1)-along(j,z))./zz(z-1);
    end
end

stress(1:56,1:84)=nan;
stress(:,2:end)=rho_constant.*amv_lat(:,2:end).*along_z;

%surface wind stress
cd /Volumes/Samsung_T5/stress_fresh
fid1=fopen('TAUX.grd','r');
fid2=fopen('TAUY.grd','r');
taux=fread(fid1,800*800,'real*4');
tauy=fread(fid2,800*800,'real*4');
fclose(fid1);fclose(fid2);
taux=reshape(taux,800,800);
tauy=reshape(tauy,800,800);
taux(taux<-10000)=nan;
tauy(tauy<-10000)=nan;
taux1=taux(r1,r2)';
tauy1=tauy(r1,r2)';

for j=1:num_size1
    taux_lat(j,1)=taux1(test3(j,1),test4(j,1),1);
    tauy_lat(j,1)=tauy1(test3(j,1),test4(j,1),1);
end
along_tau(:,1)=taux_lat(:,1).*ax_lat(:,1)+tauy_lat(:,1).*ay_lat(:,1);
stress(:,1)=along_tau(:,1);

ratio(:,:)=taup_a./stress;

for z=1:84
    along_tau_new(:,z)=stress(:,1);
end

%bottom friction stress
for i=1:56
    [along_b(i,1),num_bottom(i,1)]=bottom(i,along);%bottom velocity
    [speed_b(i,1),num_bottom1(i,1)]=bottom(i,speed);
end
btmfrc=1.3e-3;
bottom_stress=btmfrc*rho_constant*speed_b.*along_b;

for i=1:56
    stress_1(i,1:num_bottom(i,1))=stress(i,1:num_bottom(i,1));
    stress_1(i,num_bottom(i,1)+1)=bottom_stress(i,1);
    stress_1(i,num_bottom(i,1)+2:85)=stress(i,num_bottom(i,1)+1:84);
end

f=1e-4; %s-1
for z=1:84
    M_e(:,z)=(stress_1(:,1)-stress_1(:,z+1))./(f*rho_constant);
end
% the maximum stress in each column
depth_c=19;%13:90m 24:200m 19:154m
fig_scale=[39:52];
pick_stress=stress_1(fig_scale,1:depth_c)';
[M,I]=max(abs(pick_stress));

%%
cd /Volumes/Samsung_T5/paper_figures
[off_intep,depth_intep]=meshgrid(off1_dis(2:15,1),depth(1,1:depth_c));
depth_intep(1,:)=0;
figure(num_loop)
subplot(1,2,1);
pcolor(off_intep,depth_intep,stress_1(fig_scale,1:depth_c)');shading interp;
colormap(cmocean('delta'));
col1=colorbar;
caxis([-0.2 0.2])
hold on;
[ccc,cch]=contour(off_intep,depth_intep,stress_1(fig_scale,1:depth_c)',[-0.11 -0.08 -0.05 -0.02 0 0.02 0.05],...
'showtext','on','linewidth',2,'color','k');
clabel(ccc,cch,'fontsize',15,'labelspacing',200,'FontName','Arial Black','color','k');

hold on;
[ccc_1,cch_1]=contour(off_intep,depth_intep,ratio(fig_scale,1:depth_c)',[0,0],...
'showtext','on','linewidth',3,'color',[204 153 0]/255,'linestyle','-');
hold on;
[ccc_2,cch_2]=contour(off_intep,depth_intep,ratio(fig_scale,1:depth_c)',[1,1],...
'showtext','on','linewidth',3,'color',[255 99 71]/255,'linestyle','-');
hold on;
[ccc_3,cch_3]=contour(off_intep,depth_intep,ratio(fig_scale,1:depth_c)',[5,5],...
'showtext','on','linewidth',3,'color',[255 165 0]/255,'linestyle','-');

clabel(ccc_1,cch_1,'fontsize',15,'labelspacing',300,'FontName','Arial Black','color',[204 153 0]/255);
clabel(ccc_2,cch_2,'fontsize',15,'labelspacing',300,'FontName','Arial Black','color',[255 99 71]/255);
clabel(ccc_3,cch_3,'fontsize',15,'labelspacing',300,'FontName','Arial Black','color',[255 165 0]/255);

view(0,-90);
set(gca,...
    'ytick',[0 50 100 150],...
    'yticklabel',{'0','50','100','150'},...
    'FontName','Arial Black','fontsize',15,'linewidth',5,...
    'Tickdir','out');
%set(col1,'location','southoutside','position',[0.2 0.25 0.43 0.06]);
set(col1,'location','southoutside','position',[0.15 0.21 0.2 0.08]);
title(col1,'\tau(z)(kg/(m·s^2))','fontsize',20,'FontName','Arial Rounded MT Bold','color','k');
text(5,80,'no-Amur','FontName','Arial Rounded MT Bold','fontangle','italic','color','k','fontsize',20);
%print('withamur_stress(z)','-dtiff','-r600');


subplot(1,2,2);
pcolor(off_intep,depth_intep,M_e(fig_scale,1:depth_c)');shading interp;
colormap(cmocean('delta'));
col=colorbar;
caxis([-1 1])
hold on;
[ccc,cch]=contour(off_intep,depth_intep,M_e(fig_scale,1:depth_c)',[0.4 0.2 0 -0.2 -0.4 -0.5],...
'showtext','on','linewidth',2,'color','k');
clabel(ccc,cch,'fontsize',15,'labelspacing',200,'FontName','Arial Black','color','k');

view(0,-90);
set(gca,...
    'ytick',[0 50 100 150],...
    'yticklabel',{'0','50','100','150'},...
    'FontName','Arial Black','fontsize',15,'linewidth',5,...
    'Tickdir','out');
%set(col,'location','southoutside','position',[0.2 0.25 0.43 0.06]);
set(col,'location','southoutside','position',[0.59 0.21 0.2 0.08]);
title(col,'M_e(m^2/s)','fontsize',20,'FontName','Arial Rounded MT Bold','color','k');
%print('with_amur_Ekmantransport(z)','-dtiff','-r600');
set(gcf,'unit','centimeters','position',[0,0,40,14]);
%eval(['subt=suptitle(','''',name_file{num_name_file},' discharge (',name_month{num_month},' mean)','''',');']);
%
%set(subt,'FontName','Arial Rounded MT Bold','color','k','fontsize',20);
% xlabel('offshore-distance(km)','position',[-15 130],'fontname','Arial Rounded MT Bold','fontsize',20);
% ylabel('layer-depth(m)','position',[-120 -150],'fontname','Arial Rounded MT Bold','fontsize',20);
print('s3_wihout-amur','-dtiff','-r600');
num_loop=num_loop+1;

    end
end
%%
% figure(3)
% pcolor(off_intep,depth_intep,ratio(fig_scale,1:depth_c)');shading interp;
% colormap(cmocean('balance'));
% col=colorbar;
% caxis([-5 5])
% hold on;
% [ccc,cch]=contour(off_intep,depth_intep,ratio(fig_scale,1:depth_c)',[-1 0 1],...
% 'showtext','on','linewidth',2,'color','k');
% clabel(ccc,cch,'fontsize',10,'labelspacing',200,'FontName','Arial Black','color','k');
% view(0,-90);
% set(gca,...
%     'ytick',[0 50 100 150],...
%     'yticklabel',{'0','50','100','150'},...
%     'FontName','Arial Black','fontsize',15,'linewidth',5,...
%     'Tickdir','out');
% set(col,'location','southoutside','position',[0.2 0.25 0.43 0.06]);
% title(col,'\tau_p/\tau(z)','fontsize',25,'FontName','Arial Rounded MT Bold','color','k');
% %print('noriver_ratio(z)','-dtiff','-r600');

%%
function [along_b,num]=bottom(i,along_lat)
    for z=2:84
        if isnan(along_lat(i,z))
            along_b=along_lat(i,z-1);
            num=z-1;
            return
        end
    end
end


% function [ver,kkk]=transport(ver_1,ii,zz,ahv_lat)
% 
% ver=0;
% for kk=2:84
%     if ahv_lat(ii,kk)>20
%         cross_ave=(ver_1(ii,kk)+ver_1(ii,kk-1))./2;
%         trans=cross_ave*zz(kk-1);
%         ver=ver+trans;       
%     else
%         return
%     end
% end
% 
% end
% 
% 
% function index_ver=index1(ver,off_dis,int_n,hor_scale_s,hor_scale_e)
% 
% index_ver_1=0;
% for ii=hor_scale_s(int_n,1)+1:hor_scale_e(int_n,1)
%     index_ver_ave=(ver(ii,1)+ver(ii-1,1))./2.*off_dis(ii-38,1);
%     index_ver_1=index_ver_1+index_ver_ave;
% end
% index_ver=index_ver_1./sum(off_dis(hor_scale_s(int_n,1)+1-38:hor_scale_e(int_n,1)-38,1));
%    
% end